Thermal denaturation of an helicoidal DNA model 
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We study the static and dynamical properties of DNA in the vicinity of its melting transition, i.e. 
the separation of the two strands upon heating. The investigation is based on a simple mechanical 
model which includes the helicoidal geometry of the molecule and allows an exact numerical evalua- 
tion of its thermodynamical properties. Dynamical simulations of long-enough molecular segments 
allow the study of the structure factors and of the properties of the denaturated regions. Simu- 
lations of finite chains display the hallmarks of a first order transition for sufficiently long-ranged 
stacking forces although a study of the model's "universality class" strongly suggests the presence 
of an "underlying" continuous transition. 
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I. INTRODUCTION 



The study of simple dynamical models describing var- 
ious features of DNA dynamics has interested many au- 
thors for almost 20 years 0. This vivid interest arises 
of course from the biological relevance of DNA but also 
from its physical properties which can now be probed 
through single-molecule micromanipulation experiments 
like stress-induced transitions or strand separation [3| . 
This series of studies has clearly pointed out that DNA 
must be considered as a dynamical object, whose (nonlin- 
ear) distortions could play a major role in its functions. 

One feature of DNA that attracted a lot of attention 
from physicists is its thermal denaturation, i.e. the tran- 
sition from the native double-helix B-DNA to its melted 
form where the two strands spontaneously separate upon 
heating Q, because it provides an example of a one- 
dimensional phase transition. Experiments show that 
this transition is very sharp, which suggested that it 
could be first order and this led to numerous investi- 
gations, first to justify the existence of a transition in a 
one-dimensional system and second to determine its or- 
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der H H 011 H US m [13. In spite of these efforts 
the nature of the transition is still unclear and moreover, 
as discussed below, the determination of its "true" order 
from experiments may turn out to be virtually impossi- 
ble. 

In the theoretical approaches, the level of complexity 
is reduced to the minimum by taking into account only 
the (classical) motion of large subunits rather than the 
full (quantum) many-body dynamics of all the atoms. 
Clearly, an appropriate choice of the relevant degrees of 
freedom, depending on the specific problem at hand, is 
crucial. Models based on the theory of polymers use self- 
avoiding walks to describe the two strands 0, 0, ^3 ■ 
They can be very successful in studying the properties 
of the melting transition at the largest scale but, as they 
do not describe DNA at the level of the base pairs, they 
cannot be used to investigate properties that depend on 
the sequence, or probe DNA at a microscopic scale such 
as some recent single-molecule experiments. One of the 
simplest models that investigates DNA at the scale of a 
base pair is the Peyrard-Bishop (PB) model [Hill El. 
The complex double-stranded molecule is described by 
postulating some simple effective interaction among the 
bases within a pair and along the strands. The model has 
been successfully applied to analyze experiments on the 
melting of short DNA chains . Furtherm ore, it allows 
to easily include the effect of heterogeneities [I3 yielding 
a sharp staircase structure of the melting curve (number 
of open base pairs as a function of the temperature T) ^ . 
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Beyond its original motivation to explain the denatura- 
tion transition, the PB model has an intrinsic theoretical 
interest as one of the simplest one-dimensional systems 
displaying a genuine phase transition . 

The PB model has however a serious shortcoming be- 
cause it does not take into account the helicoidal struc- 
ture of DNA. In the following, we consider an extended 
version of the PB model that has been proposed to avoid 
this weakness [H Ell. The introduction of DNA geom- 
etry induces an important coupling between base pair 
opening and helical twist, largely substantiated for real 
DNA 18]. A modified version of this helicoidal model has 
also been successfully applied to describe the denatura- 
tion of the chain induced either thermally or mechani cally 
by applying an external torque to the chain ends [20j. 
Furthermore, its nonlinear excitations have been stud- 
ied: small amplitude breather-like solutions have been 
analytically determined as well as large localized 
bubbles |23| . Both types of excitations are of interest as 
they may be thermally excited as precursors for the DNA 
strand separation. A first effect of the helicoidal struc- 
ture, namely to bring close to each other bases which 
are not consecutive in the sequence, was treated in _23J, 
by introducing an interaction between these bases. This 
first approach however did not take into account the im- 
portant geometrical effect that we want to examine here, 
the coupling between opening and twist. 

The aim of this work is to give further insight on the 
melting transition of the helicoidal model, both from the 
statistical and from the dynamical point of view. After 
having recalled the model and its state variables (Sec- 
tions m and IIIII) , the first part (Sections IIVI and 0) is 
devoted to its exact thermodynamics and to a simu- 
lation study of its statistical properties. We find that 
the melting transition is extremely sharp, bearing essen- 
tially all of the hallmarks of a first-order transition, at 
all temperature sampling steps studied (down to 0.01 K). 
However, a study of the model's "universality class" , us- 
ing finite-size scaling techniques, allowing some variation 
of the relevant physical parameters, and drawing from 
analogies with a Schrodinger-like equation (Appendix A), 
strongly suggests the presence of an "underlying" con- 
tinuous transition of the Kosterlitz-Thouless type in the 
absence of nonlinear stacking. This should be contrasted 
with the exact second-order transition obtained in the 
zero-stacking limit of the PB model 9] . 

In the second part (Section IVI|) we investigate dy- 
namical structure factors that are of p articular relevance 
to both neutron and Raman [23 scattering experi- 
ments. The underlying idea is that, upon approaching 
the transition, the molecule should display precursor ef- 
fects in the form of a "mode softening", i.e. a slowing- 
down of the dynamics with appearance of a low- frequency 
component in the spectrum. Slowing effects have been, to 
some extent, observed in the structure factors of the PB 
model [T5I I thus encouraging this type of investigation. 



II. THE HELICOIDAL MODEL 

Two versions of the model have been proposed in ear- 
lier studies. The first one ^3 describes an elastic back- 
bone and fixed base-pair planes while the second con- 
siders a rigid backbone and moving base-pair planes. The 
two models display a very similar behavior with respect 
to denaturation as the potentials associated to the base- 
base interactions in a pair and along the strands are the 
same in both cases, and because both introduce a cou- 
pling between opening and twist that results from the 
helicoidal geometry. In the following, we will consider 
the fixed-planes model, sketched in Fig. 




FIG. 1: Schematic representation of the fixed-planes DNA 
helicoidal model. 

The helical structure of DNA is introduced essentially 
by the competition between a stacking interaction that 
tends to keep the base-pairs close to each other (given 
by the fixed distance h between the base planes) and the 
length £0 > h oi the backbone segment (described as an 
elastic rod of rest length ^o) that connects the attachment 
points of the bases along each strand. The ratio £o/h fixes 
the strand slant and therefore the resulting helicity of the 
structure. This helicity is accounted for by the angle of 
rotation of a base-pair with respect to the previous one, 
namely the equilibrium twist angle 9, equal to 27r/10.4 in 
B-DNA at room temperature. 

In the model, bases are described as point-like parti- 
cles of equal masses m joined by elastic rods along each 
strand. Bases lying on the same plane are coupled by 
hydrogen bonds leading to an attractive force that tends 
to maintain their equilibrium distance equal to the DNA 
diameter 2i?o. We assume that the two bases in each 
pair move symmetrically. To describe base-pair open- 
ing and helical torsion under those constraints, it suf- 
fices to introduce two degrees of freedom per base-pair: 
these are r„ and 0„, i.e. the radial and angular positions 
and of the n— th base. The total number of degrees of 
freedom is thus 2N, N being the total number of base 
pairs. The restriction imposed by the assumption of a 
symmetric motion preserves the essential feature of the 
helical structure, the coupling between torsion and open- 
ing, while it keeps the model sufficiently simple to allow 
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an exact treatment of its thermodynamics. 
We consider the Lagrangian |2q 



2 I 2 1.2 
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- I?^(exp[-a(r„-i?o)]-l 

n 
n 

- S ^(r„ - r„_i)^ exp[-&(r„ - 



2i?o)] (1) 



where = ^ + 4_Rq sin (9/2) and ln,n-i are respec- 
tively the equihbrium and the actual distance between 
the two bases n and n — 1 along a strand, 



In 



2r„_ir„ cosf 



(2) 

For later purposes, it is convenient to introduce the local 
twist angle defined as 9n = 4>n — 4>n-i- 

The first term in the Lagrangian is the kinetic energy. 
The second term is intended to describe the hydrogen 
bond interaction between the two bases in a pair. Fol- 
lowing Refs. and 0|, a simple Morse potential form 
is chosen. The quadratic term in (/„^„_i — ^o)^ represents 
the elastic energy of the backbone rods between neigh- 
boring base-pairs on each strand. Finally, the last term 
models a stacking interaction between neighboring base 
pairs. Its effect is to decrease the stiffness of the open 
parts of the chain relatively to the closed ones and to 
stabilize the latter with respect to the denaturation of a 
single base-pair. Terms of this type increase the cooper- 
ative effects close to the melting transition |23| • 

In the present paper we restrict our attention to a chain 
with free boundary conditions, which corresponds to the 
experimental situation when DNA denaturation is stud- 
ied in solution. However, the model described above can 
be easily modified to account for an external torque F 
applied at the base pairs at the two ends as it is 
done in some single-molecule experiments. 

The geometrical parameters of the model can be 
straightforwardly fixed according to the available struc- 
tural data ■ Much more delicate is instead the choice 
of the parameters b, D,S and K gauging the effective 
forces. We selected values similar to those previously 
considered for the fixed-planes case, which have been 
discussed elsewhere: the choice of the K parameter can 
be independently derived from the twist persistence 
length [23, while the choice of the other two parame- 
ters was based '29'| on a comparison with recent me- 
chanical denaturation experiments In particular, 
the important parameter D which sets the main energy 
scale, has been tuned to reproduce as closely as possible 
the experimental value of the denaturation temperature 
Td = 350 K. The full set of parameters used in the fol- 
lowing are summarized in Tabled 



Parameter 


Symbol 


Value 




Morse potential range 


a 


6.3 




Stacking interaction range 


b 


0.5 


A-i 


Morse potential depth 


D 


0.15 


eV 


Stacking interaction coupling 


s 


0.65 


eVA-2 


Inter-plane distance 


h 


3.4 


A 


Elastic coupling 


K 


0.04 


eVA-2 


Equilibrium distance 


Ro 


10 


A 


Twist angle 


e 


0.60707 


rad 


Base masses 


m 


300 


a.m.u. 



TABLE L The parameter set used throughout the paper 



For the numerics it is convenient to work in dimension- 
less units. A suitable choice is to measure lengths and 
energies in the natural units of the Morse potential, 
and D respectively, whereby time is expressed in units 
of yJmjDo?- . With the parameters of Table I, one time 
unit (t.u.) ~ 2.3 ps. 



III. THE DENATURATION TRANSITION: A 
QUALITATIVE DISCUSSION 

Before going on, it is instructive to briefly discuss the 
thermodynamic state variables of the chain as well as the 
(possible) analogies between denaturation and the more 
familiar liquid-gas transition. 

As already pointed out jl^l , the applied torque F plays 
the role of the pressure P for the liquid-gas system. Its 
conjugate variable is the degree of supercoiling a. Since 
we do not consider the curvature of the axis of the helix, 
it reduces simply to the average twist: 



N 

E 



{(On) - 9) 

N9 



(3) 



This variable thus plays the role of the volume V. Follow- 
ing this analogy we can therefore establish the following 
correspondences : 



DNA model 
T 
F 

— (7 



liquid-gas 
T 
P 
V 



The sign of a is chosen for convenience as the degree of 
supercoiling vanishes for B-DNA and is negative for a 
partially denaturated chain (see also Eq.^ below). 

Two natural scenarios may thus be expected: an 
isothermal, torque-induced denaturation at F = Fd or 
a fixed-torque, thermally-induced one at T = Tq. For 
the liquid-gas case, these two situations correspond to 
crossing of the coexistence curve in the (P, T) plane with 
an horizontal or a vertical line, leading to the transitions 
classically described by isotherms in the (P, V) plane or 
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by isobars in the (F, T) one, respectively. In both cases, 
the presence of a constant-temperature or a constant- 
pressure/torque domain is associated to the phase coex- 
istence. 



(a) 



1 






r=o 












(c) 



rc(T) 



T=Td 



FIG. 2; A sketch of the behavior of first-order transition 
curves in the three planes defined by the DNA model pa- 
rameters —<j, r and T. Td and Td represent the transition 
torque and temperature, respectively; (a) thermally-induced 
denaturation "isobar" at zero torque; (b) coexistence curve; 
(c) torque- induced denaturation isotherm at T — Td, corre- 
sponding to a zero critical torque. 

Within this analogy, the transition isotherms corre- 
spond to the curves in Fig. 3 of Ref. 20^ and repro- 
duced schematically, on the (F, — cr) plane, in Fig. 
Conversely, the thermal denaturation at constant torque 
F = are correctly described in the (— cr, T) plane 
(Fig. [2Jl). Notice that, with the convention that a neg- 
ative F corresponds to an untwisting torque, both 
and Td must increase upon increasing temperature and 
torque respectively (see Fig. I^d). However the negative 
torque cannot exceed some critical value without leading 
to an instability of the helix associated to a change of the 
sign of the helicity. In the following, we will focus on the 
thermal denaturation transition at F = 0. 

It is important to remark that the helical constraints 
included in the model roughly impose, at vanishing ex- 
ternal torque, 9n ~ for a closed chain segment, and 

« for the denatured one. This follows from the ge- 
ometry of the helix and the stiffness of the strands: since 
the distance between consecutive bases is constrained by 
the elastic rods to be approximately equal to £o, On is 
of order ^o/^n and hence very small for r„ 3> i?o- Let 
us denote by Ud the average number of open bases in a 
chain of length N . Provided that open and closed re- 
gions coexist along the helix, and that they are spatially 
well separated (this is well confirmed by simulations as 
we show below), then, to a good approximation, we have 



' N9 " N 



(4) 



The latter quantity in nothing but the average fraction of 
open base pairs p — Ud/N and the "isobar" —a{T) can 
thus be identified with the familiar denaturation curve 
p{T). In other words, the supercoiling and the fraction 
of open base pairs are equivalent order parameters. 

Before going further there is one crucial issue that 
should be addressed. One may argue that for a one- 



dimensional model like the one at hand no phase transi- 
tion of any type should be observed. However, all usual 
arguments against the existence of singularities in ther- 
modynamic potentials have been showed not to hold for 
the PB model Since the latter is in many respects 
similar to the helicoidal model, the same arguments apply 
and a genuine phase transition is not forbidden a priori. 
This is confirmed by the transfer integral approach which 
can be carried exactly (although partly numerically) for 
this simple model. 



IV. TRANSFER INTEGRAL APPROACH 
A. "Apparent" thermodynamics 

Because the model is one dimensional, a direct calcu- 
lation of the partition function can be performed by the 
transfer integral (TI) method, as it was done for the sim- 
pler PB model Igj. The calculation proceeds along the 
same lines, but it is more involved because the model 
has two degrees of freedom (r„ and 0„) per unit cell. 
It nevertheless reduces to a one-dimensional TI equation 
because the contribution introduced by the angular part 
can be diagonalizcd by a Fourier transform ,2^ ,2M ■ 
the calculation has already been presented in these ear- 
lier studies we do not discuss the method here and we 
confine our attention on its results, which point out new 
aspects of the transition that had been overlooked. 

In this subsection we restrict ourselves to results ob- 
tained via numerical solution of the TI equation. The 
accuracy of this approach is limited first by discretiza- 
tion errors in the integrations and by the need to nu- 
merically evaluate integrals over an infinite domain. As 
discussed in the next subsection, this second restriction 
can be partially lifted by a finite size scaling analysis, 
which involves a properly controlled approach to infin- 
ity. Integration methods, such as the Gauss-Legendre 
quadrature which select appropriate absissa for the eval- 
uation of the function according to the number of points 
involved in the calculation, are also useful to integrate 
over a large domain with a reasonable number of points. 
For this first study, in order to ensure that all integrals 
are evaluated with the same discretization error, we have 
computed them with a 10*'' order Bode's method [s^l 
with a fixed spatial step Sr = 0.032 A and a minimum 
value rmin = 9.7 A (due to the strong repulsion between 
bases described by the Morse potential, r cannot take val- 
ues significantly below the equilibrium length of 10 A). 
The maximum Vmax of the integration range depends on 
the number of integration points, which has been varied 
from 631 to 3601 leading to 29.9 < r„,ax < 124.9 A. The 
eigenvalues of the transfer integral operator have been 



eqm 

trix problem or by the Kellog's method Isy to get the 



obtained either by diagonalization of the equivalent ma- 
trix problem or by the 
two lowest eigenvalues. 

The eigenvalues A of the TI operator will henceforth be 
written as A = exp(— e/fc^T) where ks is the Boltzmann 
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constant. With this notation the free energy per particle 
is determined by the smallest e-eigenvalue eg and is given 
by 



(a) 



/ = —kBTln[4:TTmkBT) + eq • 



(5) 



Relevant thermodynamic quantities like the entropy and 
specific heat are then evaluated from the standard rela- 
tions 



__dl 

the mean base-pair stretching is given by 



r|0o(?')|^ dr , 



(6) 



(7) 



where (j)Q is the eigenfunction associated with eo- 

The free energy / and the mean value of the base pair 
stretching (r) for our model are shown in Figure O for 
temperatures going from 349 to 352 K with a step of 
0.02 K. Within the accuracy of the calculation, a cusp in 
/ at T]j — 350.74 K is distinctly seen. It is associated 
with a sharp jump of the entropy at the transition. A 
jump in the specific heat is also observed. Evaluating 
numerically the first and second left- and right- deriva- 
tives of the free energy, one obtains the jump in entropy 
As = 4.40 ks, or 8.75 cal/K/mol, and the jump in spe- 
cific heat Act, = 0.64 fc^. The specific heat drops from 
2.14 fcs below To to = 1.5 ks for T > Td as expected 
from equipartition because after denaturation only the 
harmonic contributions of the hamiltonian stay signifi- 
cant. 

Fig-Eb shows that, within numerical accuracy, (r) ex- 
hibits a discontinuous transition from a finite constant 
value (very close to the equilibrium value i?o) to a value 
of the order of the system size; in other words, the eigen- 
function appears to become suddenly delocalized. The 
picture of a sharp transition persists down to a temper- 
ature sampling of 0.01 K. 

Although the numerical results strongly suggest the oc- 
currence of a first order transition, caution is necessary: 
previous studies of the related PB model have shown that 
the nonlinear stacking produces an extremely sharp, first- 
order-like behavior which masks the underlying second- 
order transition as long as one stays out of a very nar- 
row domain in the immediate vicinity of the critical tem- 
perature Td (exponential crossover [s^, 113 )• A more 
complete picture of the properties of the transition will 
therefore be given in what follows. 



B. The "underlying" transition 

We first address the question of what happens in the 
absence of the nonlinear stacking, i.e. at 5 = 0. Prelimi- 
nary numerical investigations suggest a smooth behavior, 
both of the lowest eigenvalue eo and of the next-to-lowest. 




T 
(b) 



FIG. 3: The free energy per particle / (a), and the mean 
value of the base pair stretching (r) (b) of the model evalu- 
ated by the transfer integral method in the temperature range 
349 < T < 352 K with a step AT = 0.02 K. The model param- 
eters are those listed in Tabled] (Calculation with the Bode's 
method, 5r = 0.032 A, r™^i„ = 9.7 A, r^^^ = 124.9 A). 



ei, as functions of temperature; an "avoided crossing" be- 
tween them appears, with a small, but finite gap which 
has a minimum at a certain temperature. Before we pro- 
ceed to analyze the data obtained in detail, it is necessary 
to provide some background and notation. 

The order of the phase transition of the ideal system of 
unconstrained transverse spatial extent is determined by 
the critical exponent which characterizes the gap Ae = 
ei—co (X {Td^TY at temperatures below To; a value v — 
1 implies a cusp in the free energy and a discontinuous 
entropy; a value equal to 2 implies a discontinuity in the 
specific heat, i.e. a usual 2nd order transition, etc. 

The "raw" data provided by numerical solution of the 
TI equation refer to a particular transverse system size 
L = Traax determined by the imposition of an upper cut- 
off to the integration. On the other hand, near a crit- 
ical point of the infinite system, the transverse fluctu- 
ations of the order parameter also diverge. The quan- 
tity C-L = -y/ (r^) — (r)2 provides a measure of the dis- 
tance from the critical point with dimensions of length. 
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According to the finite-size scafing hypothesis |3J|, size- 
dependent properties of a system in the vicinity of the 
transition should depend solely on the ratio L/£^±, i.e. 



(8) 



where the exponent a characterizes the rounding of the 
gap (cf. below), /g(0) is a nonzero constant, and fg{x) oc 
x"^ as a; 3> 1 guarantees size-independence in the limit 
L — > cx). In the simplest cases, the positions of the min- 
ima of the gap are related to the type of divergence of ; 
according to the above scaling scenario, the temperature 
Tjn{L) where the gap minimum occurs, is such that 



e± {Td - TraiL)) ^ L 



(9) 



We have numerically solved the TI equation for a wide 
range of system sizes, using Gauss-Legendre quadratures, 
where L is defined as the largest grid point value provided 
by the Gauss-Legendre algorithm for the number N of 
grid points chosen for the calculation. Again, in order to 
ensure uniform accuracy, N grows proportionately to the 
system size, with rmax = 350 A corresponding to = 
2048 points. Fig.Qlshows the dependence of the minimal 
gap AeL(Tm), and the corresponding temperature r,„ on 
L. The gap appears to depend quadratically on 1/L (i.e. 
CT = 2 in (jSl ), to very good accuracy, with an estimated 
limit of 1.4 10~^ (with a standard deviation 0.8 10~^) as 
L oo. We note in passing that the disappearance of 
the gap between the bound state eigenvalue and the one 
belonging to the bottom of the continuum band provides 
numerical evidence for the true occurrence of an exact, 
thermodynamic phase transition [s^- The temperatures 
which correpond to the gap minima can be well fitted to 
the function 



T^{L)=Td 



02 



(10) 



with Td = 577.8A:, 02 = 0.170 and c = 0.94. This type 
of dependence of on L immediately suggests (cf. Eq. 
El) that 



a2/\r\ 



(11) 



where t = T/Td - 1. 

Fig-Elshows that data taken from a wide range of sys- 
tem sizes scale well if plotted according to Eqs. (|H)l and 
((TTIl . The numerical evidence thus strongly suggests that 
the underlying transition manifests itself as an essential 
singularity of the gap, of the Kosterlitz-Thouless (KT) 
type. Eqs. (|11|) and ||SJ| then imply that, in the limit 



L 



-2a2/|T 



(12) 



In the Appendix, it will be possible to identify the origin 
of this particular behavior as an inverse-square attractive 
interaction between the streching coordinates of succes- 
sive base pairs. 
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FIG. 4: Dependence of the minimal difference between the 
lowest eigenvalues of the transfer integral operator Aei(Tm) 
(open squares, left y-axis) and temperatures Tra (closed cir- 
cles, right y-axis) versus . 



C. Finite stacking revisited 

It is now reasonable to conjecture, by analogy with 
what happens in the PB model, that the effects of the 
nonlinear stacking interaction will depend on its range. 
For the standard parameter set of Table the ratio 
b/a = 0.079 is very small indeed. What happens at a less 
extreme regime, bja — 0.190, is shown in Fig. Scaling 
according to the Ansatz (|ll|l holds within a fairly nar- 
row range |t| < 0.05 around the dcnaturation point; note 
that it is the smaller magnitude of the nonuniversal pa- 
rameter 02 which is responsible for the narrowing of the 
asymptotic critical region. Fig. ^summarizes what hap- 
pens at 5 = 0.9 A~^, i.e. h/a = 0.143, only slightly above 
the value of Table The gap exhibits an apparent criti- 
cal exponent very close to unity down to t = 0.01; closer 
to the denaturation point, the effective slope increases 
significantly; it is reasonable to conjecture that at tem- 
peratures even closer to T^i, the asymptotic behavior will 
be dominated by the underlying essential singularity. At 
the physically relevant value of 6 = 0.5 A~^, crossover to 
the KT regime has moved below r = 10~^ and is practi- 
cally unobservable. 

The analysis of this section demonstrates that, in spite 
of the very different mathematical properties of their 
"bare" versions, both the "straight" (PB) and the he- 
licoidal DNA models, are effectively dominated by the 
stacking interaction when the latter is of sufficiently long 
range; because of it, for all practical purposes, the transi- 
tion has all the characteristics of a first order transition, 
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L/(cR„) exp(-a^/|T|) 



FIG. 5: Scaling of the difference between the two lowest 
eigenvalues of the transfer integral operator (open symbols, 
left y-scale), and the order parameter (solid symbols, right 
y-scale). The different points have been obtained by trans- 
fer integral calculations performed with 5 values of L/Rq — 
16, 35, 70, 100, 140 . The dotted lines have slopes 2 and -1, 
respectively, in accordance with the finite-size scaling hypoth- 
esis. 

including a practically infinite discontinuity of the mean 
base pair stretching, and a latent heat. Similarly, at the 
transition temperature, a very small temperature gradi- 
ent (of the order of the width of the transition region, i.e. 
AT < 0.001 K) leads to an apparent phase coexistence 
and hence to many features that one would be tempted 
to qualify as "typical of a first-order transition" as shown 
in the next sections. These properties are very reminis- 
cent of some results found on models of martensitic phase 
transitions [s^ [s^, HI] , but because we are dealing with 
a one-dimensional model that has a genuine phase tran- 
sition, the phenomenon is more remarkable here. 

Our results provide an excellent example of the dis- 
tinction between the "experimental" perspective and the 
"theoretical" one, regarding the definition of the order of 
a phase transition: although, in theory, the transition of 
the helicoidal DNA model is of infinite order (essential 
singularity), the actual temperature range over which it 
manifests its continuous character is far beyond the limits 
of either experimental or numerical observation. 

V. MOLECULAR DYNAMICS 

In this section we report the results of direct simula- 
tions of the model. They bring complementary informa- 
tions on the nature of the transition and allow us to study 
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FIG. 6: Scaling of gap (open symbols, left y-scale), and order 
parameter (solid symbols, right y-scale); model parameters 
are those of Table|l] with the exception of 6 = 1.2 A~^. The 
different points have been obtained by transfer integral cal- 
culations performed at L/Ro — 25, 35, 50. The dashed lines 
have slopes 2 and —1, respectively, in accordance with the 
finite-size scaling hypothesis. 



its dynamics as discussed in Sec. IVII As said above, we 
consider the case of thermal denaturation for F = and 
free boundary conditions. Microcanonical and canonical 
simulations were performed because they allow the ob- 
servation of the phase space from different view points. 

In the microcanonical ensemble, the Euler-Lagrange 
equations derived from were integrated directly with 
the standard fourth-order Runge-Kutta method with a 
small enough time-step (typically 0.02 t.u.) in order to 
insure that the relative energy drift is negligible (usually 
better than 10~^) on the time-scales of each run, i.e. 10^ 
to 10^ t.u. js^. Initially, all particles are set in their 
equilibrium positions (r„ — Rq, (j)n — nO) with random 
Gaussian distributed velocities (with zero average) in the 
radial direction. The variance of the distribution serves 
to fix the energy per degree of freedom e. The averag- 
ing of the quantities of interest is only started after a 
long enough transient to let the system equilibrate. Af- 
ter equilibration, the thermal energy fc^T is computed 
in the usual way as twice the average kinetic energy per 
degree of freedom. 

Constant-temperature (canonical) results were ob- 
tained through an extended Nose-Hoover method using 
a thermostat chain |40| which is specifically designed to 
constrain the total kinetic energy to fiuctuate around 
NksT, insuring at the same time the correct (canoni- 
cal) distribution of its fiuctuations. A chain of 3 ther- 
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FIG. 7: Dependence of the gap on the reduced tempera- 
ture T/Tm — 1 for h = 0.9 and system sizes L/Rq — 
12, 49, 17.5, 25, 35, 50, 70. The dashed and dotted lines have 
slopes 1 and 2, respectively. Apparent first-order behavior 
prevails for |r| > 0.01; closer to the critical point there a clear 
increase in the slope, before the onset of finite size rounding. 



mostats was employed with the first thermostat typi- 
cal frequency equal to the highest phonon frequency of 
the lattice, ujm = {a^D + 2K{Rq{1 - cose)/ lof]/mY/'^ . 
The integration of the corresponding equations of motion 
was again performed with a fourth-order Runge-Kutta 
scheme with typical time step of 0.01 t.u., and thermal- 
ization is achieved by a long-enough transient. Changes 
in temperature were performed in a sequential way upon 
heating the chain with a temperature ramp and relaxing 
afterward. 

As discussed in Sec. IIVI strictly speaking the tran- 
sition is smooth. However the temperature range of 
the crossover region to a smooth behavior is so small 
(less than AT = 0.001 K for a stacking parameter 
b = 0.5 A^^) that the numerical experiments, as well 
as actual denaturation experiments on DNA, show all 
the character of a first order transition. Therefore in this 
section we shall use the language of first order transitions 
which is the appropriate language to discuss the results. 

Measured caloric curves showing the temperature (in 
energy units) as a function of the energy per degree of 
freedom, kBT{e), for a chain of = 128 base pairs 
are reported in Fig. |S1 They distinctly show a flat part 
at ksT = 0.2012 13 corresponding to the temperature 
Td = 350 K that has been found to be the denaturation 
transition by the transfer integral calculation. In anal- 
ogy with a liquid-gas transition, the flat region occurring 
between cb = 0-20 D and ec = 0.64 D is thus identified 
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FIG. 8: Result of microcanonical simulations (open symbols): 
kinetic temperature as a function of the energy per degree of 
freedom for a molecule of A'^ = 128 base pairs. To show the 
convergence of averages, data for two transient durations are 
reported. Pluses and crosses respectively refer to canonical 
simulations with closed initial conditions and with the initial 
insertion of an artificial bubble of length 120 base pairs (tran- 
sient 10^ time units). Solid line is the transfer integral result 
(see text for details). 



as the curve of coexistence between the closed and the 
denatured phases. The slope of the two branches, fcsT/e 
should be equal to 2/c„ (the factor 2 appears because e, 
energy per degree of freedom is ^ of the energy per unit 
cell). The figure shows that this is in good agreement 
with the values given by the TI calculation, i.e. c„ ~ 2.2 
below Td and 1.5 above T^. We have also compared 
the melting entropy per particle. As = —A{df/dT), 
as obtained from the transfer integral calculation Asth., 
with that obtained from the microcanonical simulations 
Asnum. for the finite chain i.e. the ratio 2{e]j — eB)/T]j: 



As 



= 3.70 10"^eV/K, Asth. = 3.80 lO'^eV/K . 



(13) 

The two quantities are in very good agreement. 

In addition, the transition markedly displays a signa- 
ture of metastability and hysteretic effects. Indeed, the B- 
DNA branch extends well above To (up to about 500 K). 
Marked hysteretic effect upon heating are also observed 
for the thermostated chain (crosses in Fig. (S} . Either in 
microcanonical and canonical simulations, the system ap- 
pears to be spontaneously "trapped" into this metastable 
state for low enough energies (or temperatures) over the 
transition one. Direct inspection of the system configu- 
ration reveals that the chain is completely closed and we 
can refer to it as an overheated state. 

An undercooled branch exists as well below the denat- 
uration temperature. To detect it in the microcanonical 
scheme, we employed the following procedure. The ini- 
tial condition, in the B-phase, is evolved for a certain 
time after which the chain is "annealed" by multiplying 
all velocities by an assigned factor smaller than 1 (we set 
it equal to 0.8). The averages are thus computed after a 
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further transient (see again Fig. 0. Similar results can 
be obtained for the thermostatted chain by artificially 
imposing on standard initial conditions the presence of a 
denaturated bubble of given length i in the middle of the 
chain at temperature T > T]j. We will give more details 
on this procedure at the end of this section. 

Canonical and microcanonical results are therefore 
consistent, apart from some deviations at high tempera- 
tures, for T > Td, which can be expected because after 
denaturation the model becomes almost purely harmonic 
as the Morse potential linking the bases plays no role 
for base-pair distance r corresponding to the plateau of 
the potential, and the stacking contribution also tends to 
vanish. Therefore for T > Td a microcanonical equilib- 
rium cannot be achieved, unless we force it by averaging 
over thermalized initial conditions (for instance obtained 
by a Monte-Carlo procedure) or by a temporary switch 
to a canonical simulation during a run. 

The results were checked to be robust with respect to 
the transient duration as well as to the rate at which tem- 
perature is changed through the ramp. Alternative ther- 
malization schemes do not change the outcomes as well. 
For example, simulations where microcanonical runs are 
alternated to canonical ones, yield the same results (ex- 
cept a.t T > Td as mentioned above). In such a case 
the computed averages are microcanonical as the ther- 
mostatted dynamics only serves as a way to change the 
system energy. 

Obviously, a crucial issue is the dependence of the re- 
sults on the chain size. We observed that upon increas- 
ing the chain length up to iV = 256 or 512, the only 
difference with respect to the iV = 128 case is a slower 
convergence of the averages in the high and intermedi- 
ate energy region. Nonetheless, the coexistence line is 
practically reached within comparable simulation times. 
This is presumably influenced by the initial conditions 
and could be improved by a more sensible choice. 




FIG. 9: Fraction of open base pairs p and average twist {9n) 
from microcanonical simulations. The simulation parameters 
are the same as in Fig. IHl The circles were obtained with 
an initially closed chain while the squares refer to annealed 
initial conditions. 



In order to precise the nature of the transition from a 



view point close to experiments, we measured the average 
fraction of open base pairs p. This is indeed a quantity 
which is measurable by UV absorption. As in previous 
studies 10, 0, li^l , we consider a base pair to be open 
whenever the radial displacement r„ is larger than the 
inflection point of the Morse potential i.e. when r„/i?o > 
8 In 2, and we average the counting during the run. A 
simple reasoning shows that the order parameter p should 
obey a "lever rule" jij 

e = (1 - p)eB + pen, es < e < cd (14) 

thus implying that p{e) increases linearly between and 
1 along the coexistence line. In a similar way, we ex- 
pect that the average twist per base pair (6'„) — — 
(j)i)/N = 9{1 + a) should decrease linearly from a value 
close to 9 to 0. This is illustrated in Fig. ^ Notice once 
again the hysterctic effects. 

From Fig. ISjit is clear that the microcanonical ensem- 
ble has the merit of allowing to investigate the dynamics 
of the chain in the coexistence region. For illustration. 
Fig. ^1 shows a snapshot of the state of the chain for 
e/D — 0.5 were, from formula H14I) . we expect around 
40% of the base pairs to be denaturated. The fact that 
the chain opens at the sides is clearly caused by the free 
boundary conditions. Moreover this figure can give a 
hint on why we observe a transition that has all the 
features of a true first order transition (coexistence of 
phases, metastability) although the transition is actually 
second order. From a theoretical point of view, the or- 
der is determined in the thermodynamic limit, i.e. for 
an infinite system. This correspond to the identification, 
in the transfer integral method, of the free energy with 
the lowest eigenvalue eo, Eq. ||SJ). In numerical simula- 
tions, as well as in experiments, one is dealing with a 
finite system. When the thermodynamic transition is ex- 
tremely sharp (as it is the case for the stacking parameter 
6 = 0.5 A~^), the inhomogeneity caused by the free ends 
is sufficient to lead to an apparent coexistence of phases, 
i.e. a first-order-like transition, presumably because the 
boundary effects induce a perturbation (in particular on 
the average local torque) which is sufficient to change the 
local transition temperature by the very small amount 
which separates the domain of closed DNA from the do- 
main where the molecule denaturates. This is why the 
molecular dynamics simulations are useful to complete 
the transfer integral study, and to provide results that 
can be compared with experiments. 

Another interesting aspect which can be studied 
through simulation is the dynamics of opening events. 
This allows to look for analogies with the classical nucle- 
ation mechanisms that drives relaxation from metastable 
states at ordinary first-order transitions '4l'| . Fig. ^2 re- 
ports the distribution of the length of denaturated bub- 
bles for subsequent times during the same run of Fig. ^1 
There is a clear tendency for smaller bubbles to close 
(or merge) until only a few large ones remain. A sim- 
ilar measure in the overheated metastable phase shows 
instead that the size of bubbles is pretty small and de- 
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FIG. 10: Snapshot of the chain of 128 bps in the coexistence 
region e — 0.5 D. 



crease systematically in time. 

To further investigate this aspect, we performed sim- 
ulations in the canonical ensemble, starting from a ther- 
mahzed state at temperature T and artificially seeding 
a denaturated bubble of given length £ in the middle of 
the chain. To accomplish this, given the geometry of the 
model, we set 9n = in the central region and impose 
a triangular profile for the r„s designed in such a way 
that the resulting stress on the backbone springs is ap- 
proximatively zero. The flanking regions are initially at 
equilibrium and free from any additional supercoiling. 
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FIG. 11; Distributions of the lengths of the denaturated 
regions at different times in the coexistence region e = 0.5 D. 
To improve the statistics, histograms were cumulated over a 
time window of around 10^ time units. 



In a first series of simulations we checked that for 
T < Tu the B-DNA form is very robust with respect 
to such local perturbation: for instance, at T = 200 or 
300 K bubbles of length £ < 12 base pairs on a chain 
128 base-pair long tend to shrink and the system rapidly 



returns to its completely closed state. Very large bubbles 
~ 110 base pairs) may tend to close as well, although 
on longer times. However, several cases are found where 
perturbations as large as ^ = 120 base pairs are able to 
drive the system into the metastable, undercooled, denat- 
urated state. These experiment allow the observation of 
the second metastable state in the canonical scheme: thus 
completing the correspondence with the microcanonical 
results (see again data shown in Fig. |SI . 

The situation is, as expected, the opposite in the over- 
heated region To < T < 600 K. Here, the insertion of 
a short bubble suffices to destabilize the B-DNA form 
and let the system switch to its equilibrium state, i.e. 
the completely open chain. For instance, at T = 400 or 
500 K a bubble of length £ > 8 on a chain of 128 base 
pairs is generally enough. Approximated estimates seem 
to indicate that the minimal length of the bubble tend to 
decrease with temperature, as intuitively expected, but 
this behavior is not very systematic. Statistics over a 
very large number of events would be necessary to con- 
clude quantitatively. 



VI. DYNAMICAL STRUCTURE FACTORS 

One of the motivations for considering mechanical 
models is the possibility to probe microscopic and col- 
lective motion in different phases. In this section, we 
focus on dynamical correlation functions that usually re- 
flect different types of excitations. More precisely, we 
computed the radial structure factor 



Sr{q,uj) 



E 



rne 



i(qn-u}t) 



(15) 



and the angular one S^{q,uj) where ipn = fpn — nO is 
the angular displacement from the equilibrium position. 
Brackets denote an average over an ensemble of inde- 
pendent molecular dynamics trajectories (typically hun- 
dreds). All the results reported in this section are ob- 
tained in the microcanonical ensemble. 

Let us first consider the low temperature native phase. 
As shown in Fig. 1121 the spectral analysis display, as ex- 
pected, sharp lines at the frequencies of the two phonon 
branches uj±{q) that can be computed at T = in the 
harmonic approximation |27| (see the vertical lines in 
Fig. I12|l . Acoustic vibrations in the angular variables 
are only weakly coupled to the radial (optical) ones. In- 
teresting enough, the radial spectra also displays a large 
peak at a frequency lying in the phonon gap and indepen- 
dent on the wavenumber (i.e. the large peak at w ~ 0.7 
in Fig. II2I 1. Its origin can be traced back to the excita- 
tion of a localized surface mode. This is confirmed by 
direct inspection of the chain configuration. Actually, 
the mode is found to slowly decay in time due to non- 
linear interaction leading to a systematic decrease of its 
spectral component. 
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erties are still unclear although it is tempting to assign 
it to the slow dynamics of the bubble boundaries. 




FIG. 12: Structure factors Sr{q,u->) (a) and S^{q,Lu) (b) for 
A'' = 256 at very low energy, e = 0.01 D (corresponding to 
T = 18 K). The different curves correspond to wavenumbers 
q = 2.5676, 1.66896, 0.88356, 0.09812 (right to left). 



Upon increasing the energy, the optical branch gradu- 
ally shifts towards lower values of the frequency (soften- 
ing) and higher-harmonics appear. Furthermore, the res- 
onances in both the radial and angular peaks are substan- 
tially broadened due to increasing anharmonicity that en- 
hances the effective damping. More importantly, a large 
low-frequency component, a central peak, arises in the ra- 
dial structure function. The temperature dependence of 
Sriq^Lo) across the denaturation transition is illustrated 
in Fig. The three different energies correspond to 
T = 300, 357 and 535 K. The latter value is well into the 
metastable overheated region. For fixed q, the position 
of the central peak is unchanged upon increasing tem- 
perature but its width broadens. Furthermore, the uj^'^ 
behavior at low frequencies (see the inset of Fig. ll3|) sug- 
gests a Lorentzian line-shape. The origin of this central 
peak, also found in the simpler PB model, and its prop- 
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FIG. 13: Radial structure factors Sr{q,u}) for A'^ — 256 and 
different energies e = 0.17 0, 0.20 D, 0.30 D (solid, dashed, 
dot-dashed lines respectively). To reduce fluctuations, a 
smoothing of the data has been performed by averaging over 
10 consecutive channels. 



An even more sizable central component appears when 
closed and open form coexist (see Fig. I14II . This is ac- 
companied by a stronger coupling between angular and 
radial degrees of freedom, as manifested by the peaks 
at acoustic frequencies in Sr- The birth of large low- 
frequency components bears strong resemblance with 
heterophase fluctuations observed in other lattice mod- 
els with (pseudo) first-order transition characterized by 
large entropy barriers |3^ [s^, H^- In other words, the 
motion of the interface between the two phases should be 
responsible of the slow dynamics. 

To close this section, it is worth mentioning that a re- 
lated analysis of collective modes for the helicoidal model 
has beeir recently reported The analytical calcula- 
tions were performed at room-temperature and are based 
on the instantaneous normal-modes. At variance with 
our simulations, this approach describe the short-time 
dynamics (on a time scale of picoseconds) and a direct 
comparison is therefore not straightforward. 



VII. CONCLUSIONS AND DISCUSSION 

The study of a simplified model of DNA has proved to 
be extremely fruitful to unveil the basic features of the 
melting transition at the single-molecule level. From the 
theoretical point of view, dealing with a one-dimensional 
model (with two degrees of freedom per base pair) turned 
out to be particularly convenient as it allows an exact 
evaluation of the partition function. Indeed, the angular 
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FIG. 14: Structure factors Sr{q,ijj) and S^{q,uj) for A'^ = 256 
in the coexistence region e = 0.50 D. Vertical lines are the 
phonon frequencies (one from the acoustic and one from the 
optical branch for each q value) calculated at T = 0. Graphs 
have been arbitrarily shifted for clarity. 



variables can be eliminated by Fourier transform, yield- 
ing a more tractable one-variable transfer integral prob- 
lem. The latter cannot be solved analytically but the 
numerical and approximate results presented above pro- 
vide a complete insight on the nature of the transition. In 
particular, the finite-size scaling analysis of the transfer 
integral turned out to be essential to take into account 
the finiteness of the integration range. Such an analysis 
strongly suggests that the underlying transition is contin- 
uous, of the Kosterlitz-Thouless type. This behavior can 
be related to the existence of an effective attractive force 
which is directly connected to the helicoidal geometry be- 
cause it appears when the angular degree of freedom is 
integrated out. Qualitatively it can be understood as 
coming from the difficulty to disentangle the two helices. 
On the other hand, for physically relevant values of the 



parameters, the temperature range over which the con- 
tinuous aspect of the transition can be detected may be- 
come extremely narrow (less than AT = 0.001 K). For all 
practical purposes the transition appears to be perfectly 
sharp, and bears the hallmarks of a first order transition, 
in agreement with experiments. In our view, this result 
is remarkable and attracts the attention on how numer- 
ical or experimental observations on finite systems and 
with a limited resolution may dramatically differ from 
theoretical expectations. 

Molecular dynamics simulations confirm this apparent 
first-order character. They show hysteresis and metasta- 
bility as well as a coexistence region between an open and 
a closed "phase" and, by varying the energy density in 
the critical region, a gradual change of the volume frac- 
tion occupied by the two which is reminiscent of, say, the 
liquid-gas transitions. 

In addition to the very sharp transition found by the 
theoretical analysis, the finite-size effects certainly play 
a major role in the above phenomenology. The helicoidal 
model is more sensitive to these finite size effects than 
the "flat" PB model because the free ends allow a release 
of the torsional energy which appears when a segment of 
the chain opens. One can understand the crucial role of 
boundary effects if one considers that a finite closed loop 
of helicoidal DNA cannot denaturate at all because the 
two strands are entangled. 

The dynamics of the transition, as probed by the calcu- 
lation of the radial and angular structure factors, shows 
some prominent features such as the existence of a central 
peak that is presumably due to the slow motion of the 
denaturation bubbles. Moreover, the coupling between 
opening and twist introduces some additional spectral 
features that would deserve further investigations. 

Another point that should be reconsidered is the nu- 
cleation of denaturation bubbles. The phenomenology 
described at the end of Sec. is, at least qualitatively, 
very much reminiscent of the nucleation mechanisms that 
drives relaxation from metastable states at ordinary first- 
order transitions For instance, simulation in the 
overheated state suggest the existence of a "critical size" 
of the denaturation loops above which they become un- 
stable. Hence, metastability stems from the fact that 
small enough bubbles close relatively fast. Nevertheless, 
there are important differences that one should keep in 
mind. Indeed, in classical nucleation theory the key role 
is played by the surface tension term (proportional to 
the square of the droplet's radius), whereby in our one- 
dimensional case the bubble "surface" is independent of 
its length. The correspondence with the usual theory is 
probably due to the torsional energy, associated to the 
opening, which grows with the bubble size. This sug- 
gests that a "one-dimensional nucleation theory" could 
be developed for helicoidal DNA. 

The present study has focused on a DNA model that 
describes the molecule at the scale of the base pair. We 
think that it is relevant because it is the scale of the ge- 
netic code at which phenomena related to biological func- 
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tions occur. The helicoidal geometry itself is at this scale 
(or more precisely at the scale of a few tens of base pairs). 
On the other hand, there are other phenomena that enter 
in the statistics of DNA melting, and they are related to 
the behavior of the molecule at a much larger scale, on 
which the strands are regarded as flexible strings. Recent 
studies have shown that the entropy of the loops also con- 
tribute to lead to a first order transition, provided that 
self avoiding aspects betwen segments of the loop and 
between open regions and closed domains are properly 
taken into account [Tlj . Our approach is complementary 
to these studies and shows that the observed sharp melt- 
ing transition of DNA may have multiple origins. 



APPENDIX A: RESULTS OBTAINED VIA 
GRADIENT EXPANSION 

1. An approximate TI kernel 

In the absence of an external torque, the nontrivial 
part of the partition function of the model is given by 

r ^ 

(Al) 

where the potential energy consists of the three last terms 
in (Q, and the relative angle coordinate enters only via 
the second term. Introducing sum and difference coordi- 
nates r„ = (r„ -I- r„_i)/2, ^„ = r„ - Vn-i, it is possible 
to write the integral over On as 



Jo 



(A2) 



For the parameter values given in Tabled and T — 480 K, 
the dimensionless ratio A = kBT/KR^ is equal to 0.01; 
this allows us to use the leading-order low-temperature 
asymptotic expansion 



K Kf \ 4f2 



-1/2 



(A3) 



where k — 2i?o sin(6'/2) — 5.98 A, over the whole tem- 
perature range of interest. Note that due to the repulsive 
core of the Morse potential, the inequality 





r(A) 


.001 


30 


.001 


10 


.01 


10 




FIG. 15: exp(n) as a function of S/k for two different values 
of the dimensionless ratio A; the r-dependence is not visible 
for A = 0.001. 



follows. Within the approximations ljA3ll . (|A4|) . the par- 
tition function is dominated, in the thermodynamic limit, 
by the largest eigenvalue of the one-dimensional TI equa- 
tion 

/>oo 

/ dr'T{r,r')tl;^{r') ^ A^i:^{r) (A6) 

^0 



with 



T(r,r') 



K 



8f2 I 



1 



^gf2(5^)g-VM(r)/feBTg- 



8f2, 



where Vjv/(?') = f (1 — exp [— a (r — i?o)])^ and Vs — 
SS^exp [-2b {f-Ro)]. 

The form of fl (cf. Fig. [T^ and /or Eq. IA5h establishes 
that, in addition to the ranges 1/a, 1/26 of the Morse and 
stacking interactions, respectively, there is a third, much 
larger, characteristic length in the problem, k. Depend- 
ing on the strength of the various parameters, it may be 
possible to further simplify the general one-dimensional 
TI problem and elucidate the ensuing critical behavior. 
Two distinct cases will be considered below. 



r > k/2 



(A4) 



always holds. Furthermore, numerical evaluation reveals 
that the function is (i) almost independent of f, and 
(ii) weakly dependent on temperature (cf. Fig. [TK|l . In 

the following we will use the temperature-independent 
approximation 



MS') 



(A5) 



which misses the weak peak near 5 = k, but reproduces 
correctly the second moment, which is central to what 



2. Strong stacking interaction: transformation to 
an ODE 



A gradient expansion of (|A6|) involves: (i) introducing 



Mr + S)^Mr)+i^l{r)S + ^iZ{r)S^ , (A8) 

(ii) changing the variable of integration from r' {= r + 
6) to (5, and (iii) performing the Gaussian integrals over 
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5. Noting that (a) the combined effect of the stacking 
interaction and the Gaussian approximation (|A5I) can be 
described in terms of the quantity 



S 



-2b{r-Ro) I J_ 



(A9) 



( note that can be interpreted as an effective nearest 
neighbor harmonic spring constant ) and that (b) f « r 
to second order in S everywhere in l|A7(l . one obtains 



1 - 



16S'/i^r^ 



knT 



+ 4^"^^ ^ e-''(^-^^)V^. (AlO) 



where Ui{r) = Vnir) + VL(r) + Vsir) and = 
{TikBT/y/WsyQ/ne-^^-, here, VlW = -kBT{K/rf/8 
is a long-range attraction which comes from in expo- 
nentiating the term in the first parentheses of ljA7p . and 
VB{r) — k bT \n[^{r) / ^{oo)] is a thermally generated bar- 
rier analogous to the one described in 8] in the context of 
the one-dimensional DNA model with stacking. Expand- 
ing the exponential in the r.h.s. of (|A10p and rearranging 
terms, one obtains 



{ksTf 



+ (Vm + Vb + VI) = e^ip^ 



where 



Vlir)^VL{r) 
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(A12) 



is attractive everywhere. 

Eq. (|A11|I is a key result. It can be trivially cast 
in standard Sturm-Liouville form with a density func- 
tion proportional to /i^; the r dependence of fi is crucial 
for obtaining quantitatively sensible results; the simul- 
taneous presence of three terms in the potential energy 
prevents us from solving (jAllll exactly. A numerical so- 
lution j44j for the parameter values of Table Q reveals a 
behavior very similar to the full TI solution of Scction llVl 
According to Fig. ^| the two lowest eigenvalues exhibit 
an almost perfect intersection at a temperature sampling 
AT = 0.1 K. In addition, the differential equation turns 
out to be an excellent approximation to the original TI. 
Thus, the estimated To = 370 K is only a few percent 
higher than the value obtained within the TI; other crit- 
ical thermodynamic quantities of interest demonstrate 
comparable, or better agreement, e.g. the transition en- 
thalpy AH = TdAS = 0.129 eV (cf. 0.133 eV from 
TI), or the jump in the specific heat, O.Sfcs (cf. Q.lks 
from TI). The apparent first-order transition has its ori- 
gin in the fact that the thermally generated barrier has 
a substantially longer range than the Morse potential. 
The analysis of Section ^\ suggests that crossover to a 
continuous transition eventually occurs; however, for the 
values of the parameters relevant to DNA denaturation, 
observing this exponential crossover would require a tem- 
perature resolution of better than 1 mK. This estimate 
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FIG. 16: The numerically determined two lowest eigenvalues 
of llAllL expressed in units of D. The onset shows details of 
the gap in the region of the transition; no signs of rounding 
are apparent at a sampling of AT = 0.1 K. 



can be made by studying crossover phenomena with ex- 
actly solvable "toy models" 's^l of the denaturation tran- 
sition of the linear PB variety, where the zero stacking 
limit is known to yield a second-order " underlying" tran- 
sition. In the case of IjAlip . the presence of an attractive 
inverse square interaction raises the possibility of more 
complex behavior, i.e. confinement at all temperatures 
or crossover to another transition at much higher tem- 
peratures. This is discussed below. 



The 5 = case 



In the limit S 
governed by (|A5|) 



0, the decay of the kernel HA7|I 
This reduces HA11|I to 



-kBT'^^jj^ 



where 



-kuT 



(A13) 



(A14) 



In the absence of stacking, the system is subject to the 
Morse potential and the long-range attraction ljA14p : the 
point to note is that the attractive force is linearly de- 
pendent on the temperature, just as the coefhcient of 
the 2nd derivative in (|A13|) : consequently, if D = 0, the 
system will either have a bound state or not, accord- 
ing to the value of the coefficient in the denominator of 
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(jll4|l . The value 16 is marginal; if the interaction had 
been stronger, one would have confinement at all temper- 
atures; the Morse potential, being of short range, could 
not change that; in other words, one would obtain a near- 
transition at a temperature controlled by the Morse po- 
tential, but then the long-range attraction would prevent 
dissociation at all temperatures. We have verified this by 
numerically solving (jA13|l . For weaker attractions (value 
of the coefficient 16 or higher in these units), numerical 
work suggests that the transition becomes higher than 
second order; however, numerical accuracy is not suffi- 
cient to determine the detailed behavior. It is possible to 
guess what happens by substitutng the Morse potential 
by a narrow well, i.e. the total potential in l|A13|l being 
equal to —D for i?o < r < Rq + 1/a and equal (|A14|I 
for larger r; this case is exactly solvable and shows that 



although the shift in the value of the critical point is less 
than 1 % , the nature of the transition is radically trans- 
formed: the vanishing of the lowest eigenvalue is now of 
the Kosterlitz-Thouless type 
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